Delocalization of wave packets in disordered nonlinear chains 
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We consider the spatiotemporal evolution of a wave packet in disordered nonlinear Schrodinger and 
anharmonic oscillator chains. In the absence of nonlinearity all eigenstates are spatially localized 
with an upper bound on the localization length (Anderson localization). Nonlinear terms in the 
equations of motion destroy Anderson localization due to nonintegrability and deterministic chaos. 
At least a finite part of an initially localized wave packet wiU subdiffusively spread without limits. 
We analyze the details of this spreading process. We compare the evolution of single site, single mode 
and general finite size excitations, and study the statistics of detrapping times. We investigate the 
properties of mode-mode resonances, which are responsible for the incoherent delocalization process. 

PACS numbers: 05.45.-a, 05.60.Cd, 63.20.Pw 



I. INTRODUCTION 

The normal modes (NMs) of a d = 1-dimensional lin- 
ear system with uncorrelated random potential are spa- 
tially localized (Anderson localization). Therefore any 
wave packet, which is initially localized, remains local- 
ized for all time [l[ . Note that NMs correspond to single 
particle eigenstates of related quantum systems. 

When nonlinearities are added, NMs interact with each 
other 12] . Recently, experiments were performed on light 
propagation in spatially random nonlinear optical media 
[3, 01 and on Bose-Einstein condensate expansions in ran- 
dom optical potentials ^5}], which serve as realizations of 
such cases. 

Numerical studies of wave packet propagation in sev- 
eral models showed that the second moment of the 
norm/energy distribution grows subdiffusively in time as 
with a « 1/3 for d = 1. Reports on partial 
localization were published as well • 

In a recent letter the mechanisms of spreading and lo- 
calization were studied for d = 1, with initial excitations 
being localized on a single site [iTi. A theoretical ex- 
planation of the exponent a ~ 1/3 was obtained, con- 
sistently assuming that the internal dynamics of a wave 
packet is chaotic, leading to a partial dephasing of the 
NMs. The argumentation was based on the possibility 
of a pair of wave packet modes being able to resonantly 
interact with each other. Among other results, the case 
of weak nonlinearity showed that wave packets localize 
according to the linear dynamics on long but finite time 
scales, with subsequent detrapping. In the present work, 
we extend this study to single mode excitations, and more 
general excitations of width L. We study the details 
of the detrapping process, and measure the statistical 
properties of detrapping times. We study the particular- 
ities of resonant interaction between modes, mediated by 
the nonlinearity. We give details on the used integration 
schemes, and perform extensive tests which demonstrate 
that the observed effects are not affected by roundoff er- 
rors. We argue that the spreading is inherently induced 



by the nonintegrability of the system. 

II. MODELS 

We study two models of one-dimensional lattices. 

A. Nonlinear Schrodinger lattice 

The Hamiltonian of the disordered discrete nonlinear 
Schrodinger equation (DNLS) 



{4>i+i*; + *;+,*,) (1) 



with complex variables ipi, lattice site indices I and non- 
linearity strength /? > 0. The random on-site energies e/ 
are chosen uniformly from the interval [— ^] , with W 
denoting the disorder strength. The equations of motion 
are generated by ipi — dHo/diiipf): 



(2) 



Eqs. ([2]) conserve the energy ((!]) and the norm S = 
J2i We note that varying the norm of an initial 

wave packet is strictly equivalent to varying /3, there- 
fore we choose S — 1. Eqs. ^ and ^ are derived 
e. g. when describing two-body interactions in ultracold 
atomic gases on an optical lattice within a mean field ap- 
proximation [l2|, but also when describing the propaga- 
tion of light through networks of coupled optical waveg- 
uides in Kerr media 13]. 

For (3 = Eq. ((!]) with tpi = Ai exp{—i\t) is reduced 
to the linear eigenvalue problem 



XAi = eiAi ~ Ai^i - A 



(3) 



The normalized eigenvectors A^^i (X);^^/ = 1) ^-^^ the 
NMs, and the eigenvalues are the frequencies of the 
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NMs. The width of the eigenfrequency spectrum of 
® is Ac = + 4 with g [-2 - f , 2 + f ] . 

The asymptotic spatial decay of an eigenvector is given 
by A^^i ~ e-'/«(^") where_£(Ai,) < ^(0) « 
is the localization length [iJl- The NM participation 
number — l/X^z^^i characterizes the spatial ex- 
tend (localization volume) of the NM. It is distributed 
around the mean value ~ 3.6^(Ay) with variance 
« 1.3^(Aj,) The average spacing of eigenvalues of 

NMs within the range of a localization volume is there- 
fore AXd w Ad/pU « AnW^/SeO. The two scales 
AXd < Ad determine the packet evolution details in 
the presence of nonlinearity. 

The equations of motion of ([1]) in normal mode space 
read 

(4) 

1/1 ,1/2,1/3 

with the overlap integral 

Iu,Vx,U2,V3 — ^ Ajj^lAu-^^^lAy^^lAn^^l . (5) 

I 

The variables (f>ij determine the complex time-dependent 
amplitudes of the NMs. 

The frequency shift of a single site oscillator induced by 
the nonlinearity \s 5i = (3\ipi\'^ . If instead a single mode 
is excited, its frequency shift is given by 5i, — (3\(j)u\^ /p^. 

B. Anharmonic oscillator lattice 

The Hamiltonian of the quartic Klein-Gordon lattice 
(KG) 

= E y + + \< + ^(-'+1 - (6) 

where ui and pi are respectively the generalized co- 
ordinates and momenta, and ei are chosen uniformly 
from the interval [^,f]- The equations of motion are 
ill = —BT-Lk I dui and yield 

ui = -iiui -uf + ^{ui+i + ui^i - 2ui) . (7) 

Equations (O conserve the energy They serve e.g. 
as simple models for the dissipationless dynamics of an- 
harmonic optical lattice vibrations in molecular crystals 
[l^ . The energy of an initial state E >0 serves as a con- 
trol parameter of nonlinearity similar to (3 for the DNLS 
case. 

The coefficient 1/(2M^) in ^ was chosen so that the 
linear parts of Hamiltonians ([T]) and ^ would correspond 
to the same eigenvalue problem. In practice, for E — > 
(or by neglecting the nonlinear term wf/4) model ^ 
with ui = Ai exp{iujt) is reduced to the linear eigenvalue 
problem © with A = Wuj^ - W ~ 2 and ei = W{ei - 1). 



The width of the squared frequency spectrum is Ak = 
1 + ^ with ujI e [5, i + w] • ^^^^ = I^A^. As 

in the case of DNLS, W determines the disorder strength. 

The spatial properties of the NMs are identical with 
those of ([3]). In addition to the scale A^, the average 
spacing of squared eigenfrequencies of NMs within the 
range of a localization volume is Aw^ = Ak/pU- The 
two scales Alu"^ < Ak determine the packet evolution 
details in the presence of nonlinearity. 

The squared frequency shift of a single site oscillator in- 
duced by the nonlinearity is Si ~ (3£'/)/(2ei), where Ei is 
the energy of the oscillator. If instead a single mode is ex- 
cited, its frequency shift is given by 6„ ~ {3E,y) / {2p^ujl) 
with El, being the energy of the mode. 

For small amplitudes the equations of motion of the 
KG chain can be approximately mapped onto a corre- 
sponding DNLS model (T^j. In our notation, the map- 
ping takes the following form. For the KG model with 
given parameters W and E, the corresponding DNLS 
model ([1]) with norm S = 1, has a nonlinearity param- 
eter (3 ~ 3WE. The norm density of the DNLS model 
corresponds to the normalized energy density of the KG 
model. 



C. Computational methods 

We will present results on long time numerical simula- 
tions. We therefore first discuss the methods and partic- 
ularities of our computations. For both models, we used 
symplectic integrators. These integration schemes re- 
place the original Hamiltonian by a slightly different one, 
which is integrated exactly. The smaller the time steps, 
the closer both Hamiltonians. Therefore, the computed 
energy (or norm) of the original Hamiltonian function 
will fluctuate in time, but not grow. The fluctuations are 
bounded, and are due to the fact, that the actual Hamil- 
tonian which is integrated, has slightly different energy. 

Another possible source of errors is the roundoff pro- 
cedure of the actual processor, when performing oper- 
ations with numbers. Sometimes it is referred to as 
'computational noise' although it is exactly the oppo- 
site, i. e. purely deterministic and reproducible. We will 
discuss the influence of roundoff errors on our results in 
section IhTfI 

The KG chain was integrated with the help of a sym- 
plectic integrator of order 0(t'^) with respect to the inte- 
gration time step r, namely the SABA2 integrator with 
corrector (SABA2C), introduced in A brief presen- 
tation of the integration scheme, as well as its imple- 
mentation for the particular case of the KG lattice ([6]) is 
given in Appendix [XI The SABA2C integration scheme 
proved to be very efficient for long integrations (e. g. up 
to 10^*^ time units) of lattices having typically N — 1000 
sites (see for example the right plots of Fig. [2]), since it 
kept the required computational time to feasible levels, 
preserving at the same time quite well the energy of the 
system. For example, an integration time step r = 0.2 
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usually kept the relative error of the energy smaller than 
10-*. 

The DNLS chain was integrated with the help of the 
SBAB2 integrator (see Appendix [X)) . which introduces 
an error in energy conservation of the order ©(r^). The 
number of sites used in our computations varied from 
N = 500 to = 2000, in order to exclude finite size 
effects in the evolution of the wave packets. For r = 0.1 
the relative error of energy was usually kept smaller than 
10"'^. It is worth mentioning that, although the SBAB2 
integrator and the commonly used leap-frog integrator 
introduce errors of the same order, the SBAB2 scheme 
exhibits a better performance since it requires less CPU 
time, keeping at the same time the relative energy error 
to smaller values than the leap-frog scheme. 

We order the NMs in space by increasing value of the 
center-of-norm coordinate Xi, = J2i / • We analyze 
normalized distributions z^, > using the second moment 
™2 = J^ui'^ ~ ^Y'^v, which quantifies the wave packet's 
degree of spreading and the participation number P — 
1/ Si/ -2^5 which measures the number of the strongest 
excited sites in z^. Here V = ^^vz^. For DNLS we 
follow norm density distributions Zi, = \4>u\'^ / l^pP- 
For KG we follow normalized energy density distributions 
z^ EE E^l Ef, with = Al/2 + LjlAl/2, where is 
the amplitude of the i^th NM and ivl ^ 1 + (A,^ + 2)/W. 

III. WAVE PACKET EVOLUTION 

Below we will mainly use the DNLS case for theoretical 
considerations, and also discuss crucial points to be taken 
into account, when considering the KG case. We will 
present numerical results for both models. 

We first consider a wave packet at i = which is 
compact either in real space, or in normal mode space. 
Compactness in real space implies a single site excitation 
"0; — ^i,io with the choice eig = for the DNLS model. 
For the KG model we set pi — 0, ui — cSi^ig with iig = 1 
and c being a constant which defines the initial energy 
E. Compactness in normal mode space instead implies 
a single mode excitation 0^ — S^^^^ with X^^ w for the 
DNLS model, while in the case of the KG system we have 
=cS^,^,, = 0, withw^^ « l + (2/W), i. e. ujI^ is lo- 
cated in the middle of the frequency spectrum. Again the 
constant c defines the initial energy of the wave packet. 
We will later also consider finite size initial distributions 
of width L. 



A. Expected regimes 

Let us consider a single site initial excitation with a 
corresponding nonlinear frequency shift 61. We compare 
this frequency shift with the two scales set by the linear 
equations: the average spacing AA (which corresponds 
to AAd for DNLS and to Aw^ for KG) and the spectrum 
width A (with A denoting A^i for DNLS and Ak for 




FIG. 1: (color online) Schematic representations of the three 
different regimes of spreading for the DNLS (left graph) and 
the KG model (right graph), in the parameter space of dis- 
order strength W and of the nonlinear frequency shift S at 
initial time t = 0. For each regime the dependence of logm2 
(blue solid curves) and of log P (red dashed curves) versus 
\ogt are shown schematically (see section UlI CI for details). 



KG). We expect three qualitatively different dynamical 
regimes: I) 5i < AA; II) AA < ^; < A; III) A < Si. In 
case I the local frequency shift is less than the average 
spacing between interacting modes, therefore no initial 
resonance overlap of them is expected, and the dynamics 
may (at least for long times) evolve as in the linear case 
(/3 = for DNLS and ^ for KG). In case II reso- 
nance overlap may happen immediately, and the packet 
should evolve differently. For case III the frequency shift 
exceeds the spectrum width, therefore some renormalizcd 
frequencies of NMs (or sites) may be tuned out of reso- 
nance with the NM spectrum, leading to selftrapping. 
The above definitions are highly qualitative, since local- 
ized initial conditions are subject to strong fluctuations. 

If we instead consider a single mode initial excitation, 
we have to replace Si by Si, in the above argumenta- 
tion. For both the DNLS and the KG model, it fol- 
lows Si ^ PvSv The mean NM participation number 
(the localization volume) > 1 depends on the disorder 
strength W . 

If an initial excitation of the DNLS model is charac- 
terized by some exponentially localized (not necessarily 
compact) distribution ipi with 5 = 1, the nonlinear fre- 
quency shift may be roughly estimated as S ~ 
where the maximum norm density IV'P = sup; IV'iP- The 
left graph of Fig[T]shows the location of the three different 
regimes in the plane of the control parameters, i. e. the 
frequency shift S and the disorder strength W . Note that 

and the intermediate regime 
II disappears around VF « 20, where the participation 
number of a NM becomes of the order of one, and the 
NMs become almost single site solutions. Similarly, for 
the KG model we have the estimation S ^ E and the cor- 
responding parameter space of the three different regimes 
is shown in the right graph of FigU] 



B. The selftrapping theorem 

Regime III is also captured by a theorem presented in 
p^ . which proves, that for /3 > A (for the DNLS case) 
the single site excitation can not uniformly spread over 
the entire (infinite) lattice. Indeed, with the notations 



n 



NL 



/3 



p: 



(8) 
(9) 

(10) 



where Pr is the participation number in real space, the 
single site excitation at time t = yields 



TiLit = 0) = , HNLit = 0) 



(11) 



Due to norm conservation = 1 at all times, the har- 
monic energy part Ti^ is bounded from above and below 



0: 



(12) 



Due to energy conservation, for all times the anharmonic 
energy part Ti.NL can therefore not become smaller than 



B W 



(13) 



It follows with (|10p . that the participation number is 
bounded from above by a finite number, which diverges 
for /3 = A: 



Pr{t) < 



(3- A 



if /3 > A . 



(14) 



Moreover, since Pr ^ =J2i IV'/ 1'' < sup; \ipi\'^ [K 
elude that 

sup|^,p(^)>^- 



(15) 



Therefore, at least a part of the wave packet will not 
spread, and stay localized, although the theorem does not 
prove that the location of that inhomogeneity is constant 
in time. The norm of the part of the wave packet, which 
can spread uniformly over the entire system, is bounded 
from above by 5*00 < A//3. 



C. Numerical results 

We first show results for single site excitations [III ■ We 
systematically studied the evolution of wave packets for 
lattices ([T]) and The scenario described in section 
nil Al was observed very clearly. Representative examples 
are shown in FigO Regime III yields selftrapping (see 



10' 



10 



10 



10 

10^ 



10 




I L 



I 1 1 iiiiiii p 




I L 




I L 




I L 



mniii I 



lo' lo' lo' lo' lo' lo' lo' io'° 
t t 



FIG. 2: (color online) Single site excitations. m2 and P ver- 
sus time in log-log plots. Left plots: DNLS with W = 4, (3 = 
0,0.1,1,4.5 [(o), orange; (b), blue; (g) green; (r) red]. Right 
plots: KG with W = 4 and initial energy E = 0.05,0.4, 1.5 
[(b) blue; (g) green; (r) red]. The orange curves (o) corre- 
spond to the solution of the linear equations of motion, where 
the term uf in (O was absent. The disorder realization is 
kept unchanged for each of the models. Dashed straight lines 
guide the eye for exponents 1/3 (m2) and 1/6 (P) respec- 
tively. Insets: the compactness index ^ as a function of time 
in linear-log plots for /3 = 1 (DNLS) and E = 0.4 (KG). 



also Figs. 1, 3 in [T3|), therefore P does not grow signifi- 
cantly, while the second moment m2 ~ t"' with awl/3 
(red curves). Thus a part of the excitation stays highly 
localized [13], while another part delocalizes. Regime II 
yields subdifi'usive spreading with m2 ~ and P ~ 
0, Hi (green curves). Regime I shows Anderson localiza- 
tion up to some time which increases with decreasing 
nonlinearity. For t < Td both m2 and P are not chang- 
ing. However for t > a. detrapping takes place, and the 
packet starts to grow with characteristics as in regime II 
(blue curves). The simulation of the equations of mo- 
tion in the absence of nonlinear terms (orange curves), 
demonstrates the appearance of Anderson localization. 

The second moment m2 is sensitive to the spreading 
distance of the tails of a distribution, while the participa- 
tion number P is a measure of the inhomogeneity of the 
distribution, being insensitive to any spatial correlations. 
Thus, P and m2 can be used to quantify the sparseness of 
a wave packet. To this end, we introduce as a measure of 
the compactness of a wave packet the compactness index 



c 



p2 
7712 



(16) 



Let us consider a wave packet of K sites {K ^1). In 
the case where all the K sites are equally excited the com- 
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FIG. 3: Norm density distributions in the NM space at time 
t = 10* for the initial excitations in the regime II of the 
DNLS model shown in the left plots of Figs. [2] and [5] Left 
plots: single site excitation for W — 4 and P — 1. Right 
plots: single mode excitation for W — 4 and /3 = 5. \(j>v\'^ 
is plotted in linear (logarithmic) scale in the upper (lower) 
plots. The maximal mean value of the localization volume 
of the NMs p ~ 22 (shown schematically in the lower plots) 
is much smaller than the length over which the wave packets 
have spread. 



pactness index is given by C = 12. In the case of a sym- 
metric wave packet formed by a sequence of an excited 
site followed by a nonexcited one, where all the K/2 ex- 
cited sites have the same amplitude, C = 3. Distributions 
with larger gaps between the equally excited isolated sites 
attain a compactness index C, < 3. For the extreme case 
of a sparse wave packet formed by two equally excited 
sites located at the two edges of the packet, i. e. when 
only sites 1 and K [K ^ 1) are excited to an amplitude 
1/2, the compactness index is C = 16/ K^. So, smaller 
values of C correspond to more sparse wave packets. 

We expect that ^ in regime I will remain constant for 
t < Td and will behave as in the case of regime II for 
latter times. In regime II ( would either be constant or 
decay in time, while in regime III it should decay since P 
remains practically constant. The time evolution of C for 
excitations in regime II is shown in the insets of Fig. [2l 
As one can see the compactness index oscillates around 
some constant nonzero value both for the DNLS and the 
KG models. This means that the wave packet spreads 
but does not become more sparse. For the particular 
cases of Fig. [5] the compactness index attains the values 
C = 3.5 for the DNLS model at t = 10^ and C = 1-7 
for the KG chain aX t = 10^". The corresponding wave 
packet of the DNLS model is shown in the left plots of 
Fig.H 

Partial nonlinear localization in regime III is explained 
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FIG. 4: (color online) Single site excitations for the same 
disorder realization of the KG model, and P versus time 
in log-log plots. Left panels: plots for W — 4 and initial 
energy E = 3.225,4,10 [(bl) black; (r) red; (g) green]. Right 
panels: Plots for E = 0.05 and W = 6,7 [(bl) black; (r) red]. 



by selftrapping [T^j . It is due to tuning frequencies of ex- 
citations out of resonance with the NM spectrum, takes 
place irrespective of the presence of disorder and is re- 
lated to the presence of exact i-periodic spatially local- 
ized states (also coined discrete breathers) for ordered 
[Tot and disordered systems [l^j (in the latter case also 
t-quasiperiodic states exist). These exact solutions act 
as trapping centers. 

Note that for large nonlinearities (/? ^ A for DNLS 
or large energy values E of the KG model) almost the 
whole excitation is selftrapped. This behavior can be 
seen in the left plots of Fig. [H where the time evolution 
of 7712 and P for different values of the energy E of the 
KG chain is shown. The value of W is kept to W = 4 
as in the cases presented in the right plots of Fig. [2l 
As the energy increases the portion of the wave packet 
that stays selftrapped increases with respect to the part 
that diffuses. Thus, we observe a change in the evolution 
of m2 from subdiffusive increase to practical constancy. 
On the other hand, P is not affected as it continues to 
fluctuate around some constant value. 

Anderson localization on finite times in regime I is ob- 
served on potentially large time scales and as in III, 
regular states act as trapping centers [20]. For t > Td, 
the wave packet trajectory finally departs away from 
the vicinity of regular orbits, with subsequent spread- 
ing. Increasing the value of W results to small localiza- 
tion lengths of NMs and thus, Anderson localization will 
persist for extremely long time intervals. Since our nu- 
merical computations are limited in time, we are not able 
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FIG. 5: (color online) Single mode excitations. m2 and P 
versus time in log-log plots. Left plots: DNLS with W = 4, 
/3 = 0, 0.6, 5, 30 [(o) orange; (b) blue; (g) green; (r) red]. Right 
plots: KG with W ^ 4 and initial energy E = 0.17, 1.1, 13.4 
[(b) blue; (g) green; (r) red]. The orange curves (o) corre- 
spond to the solution of the linear equations of motion, where 
the term uf in ((Tj was absent. The disorder realization is 
kept unchanged for each of the models. Dashed straight lines 
guide the eye for exponents 1/3 {m2) and 1/6 (P) respec- 
tively. Insets: the compactness index as a function of time 
in linear-log plots for /3 = 5 (DNLS) and E = 1.1 (KG). 



to observe the detrapping phase of the evolution when W 
increases significantly. This behavior can be seen in the 
right plots of Fig. |4] where we consider initial single site 
excitations which, for W = 4 (see right plots of Fig. ^ 
belong to regime I. In these plots we observe a direct tran- 
sition from regime I to practical constancy of m2 and P 
as W increases, at least up to the final integration time 
used. 

For single mode excitations we find a similar out- 
come, but with rescaled critical values for the nonlinear- 
ity strength which separate the different regimes. Exam- 
ples of the three different regimes are shown in FigO As 
in the case of single site excitations presented in Fig. [21 
the compactness index ( plotted in the insets if Fig. [5] 
remains practically constant for excitations in regime II, 
attaining the values C = 1-5 at t = 10® for the DNLS 
model and ( — 3.3 a,t t — 10^ for the KG chain. The 
final norm density distribution for the DNLS model is 
plotted in the right plots of Fig. [3l The average value ( 
of the compactness index over 20 realizations a.t t = 10® 
for the DNLS model with W = 4 and (3 — 5 was found 
to be C = 2.95 ±0.39. 



FIG. 6: (color online) Single site excitations. m2 (in arbitrary 
units) versus time in log-log plots in regime II and different 
values of W. Lower set of curves: plain integration (without 
dephasing); upper set of curves: integration with dephasing of 
NMs (see section |IVA[I . Dashed straight lines with exponents 
1/3 (no dephasing) and 1/2 (dephasing) guide the eye. Left 
plot: DNLS, W = 4, /3 = 3 (blue); W ^ 7, P = 4 (green); 

= 10, /3 = 6 (red). Right plot: KG, W ^ 10, E = 0.25 
(blue) , W = 7, E ^ 0.3 (red) , W = 4, E = 0.4 (green). 
The curves are shifted vertically in order to give maximum 
overlap within each group. 



D. Spreading 

The subdiffusive spreading takes place in regime I for 
t > Td, in regime II, and for a part of the wave packet 
also in regime III. For single site excitations the expo- 
nent a does not appear to depend on (3 in the case of 
the DNLS model or on the value of E in the case of KG. 
In Figl6] we show results for m2 (t) in regime II for dif- 
ferent values of the disorder strength W. Again we find 
no visible dependence of the exponent a on W . There- 
fore the subdiffusive spreading is rather universal and the 
parameters /? (or E) and W are only affecting the pref- 
actor. Excluding selftrapping, any nonzero nonlinearity 
will completely delocalize the wave packet and destroy 
Anderson localization. We performed fittings by analyz- 
ing 20 runs in regime II with different disorder realiza- 
tions. For each realization we fitted the exponent a, and 
then averaged over all computational measurements. We 
find a = 0.33 ± 0.02 for DNLS, and a = 0.33 ± 0.05 for 
KG. Therefore, the predicted universal exponent a = 1/3 
[m appears to explain the data. 

On the other hand, in the case of single mode excita- 
tions the numerically computed values of the exponent 
a seem to be slightly larger than a = 1/3, as can be 
also seen from the results of Fig. O In particular, m2 
in regimes II and III of the DNLS model and in regime 
III of the KG model increases slightly faster than oc t-^^^, 
which is represented by the dashed lines in the upper 
plots of Fig. [51 In addition, the value of the exponent 
seems to slightly vary with respect to the nonlinearity 
parameter f3 for DNLS and E for KG. The reason of the 
slightly different behavior between single site and single 
mode excitations is still an open issue. 
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FIG. 7: (color online) Evolution of 7712 versus time in log-log 
plots. Single site excitations in the intermediate regime II for 
the DNLS (left plot) and the KG model (right plot) corre- 
spond to black curves (bl). The wave packets after td = 10^, 
lO"*, 10^, 10*^ time units (t. u.) [(r) red; (g) green; (b) blue; 
(p) purple] are registered and relaunched as initial distribu- 
tions (colored curves). The dashed straight lines correspond 
to functions oc t^^^. 



E. Detrapping 

In the intermediate regime II the wave packet starts to 
spread almost from scratch. We do not observe any satu- 
ration and crossover into localization on later times. Let 
us assume that the wave packet spreads without Umita- 
tions. The initial nonlinear frequency shift Si was larger 
than the average level spacing in a localization volume 
AA. However, Si will become smaller than AA at some 
later time, since sup; (sup; Ei for KG) decreases in 
time as the wave packet spreads. Therefore, there will be 
a large but finite time td, at which we cross over from the 
intermediate regime II into the weak nonlinearity regime 
I. The arresting of the wave packet up to a time in the 
weak nonlinearity regime I can be explained by a corre- 
spondingly large spreading time scale r^. For t < Td no 
spreading is observed when monitoring the second mo- 
ment m2, with subsequent spreading observed on larger 
time scales t > Td. 

We test the above conclusions by the following simple 
scheme. We start a single site excitation in the interme- 
diate regime II, measure the distribution at some time td, 
and relaunch the distribution as an initial condition at 
time t = 0. The results are shown in FiglT] We find that 
the relaunched runs yield a second moment m2 which 
appears to be constant up to the time Td ~ td with a 
subsequent spreading, similar to the previously obtained 
detrapping in regime I. 

For a specific value of the nonlinearity (3 of the DNLS 
model let each NM in the packet after some spreading to 
have norm ~ n <C 1 with n denoting the average 

norm density of the excited NMs (in the case of the KG 
model n corresponds to the average energy density of the 
excited NMs). The packet size is then 1/n ^ p, with 
p = maxj^ p^, and the second moment TO2 ^ l/n^. Let 
us assume that the second moment grows as 7712 ~ t-^^^. 
Let us also assume, that at any time the spreading is due 




FIG. 8: (color online) Nonlocal excitations of the KG chain 
corresponding to initial homogeneous distributions of energy 
E = 0.4 over L neighboring sites, (a) m.2 and (b) P versus 
time in log-log plots for L =1, 9, 19, 29 and 39 sites [(bl) 
black; (r) red; (g) green; (b) blue; (p) purple], (c) Fitting 
of the time evolution of 7712 for L = 19 with a curve of the 
form (HHI for M = 3.25, Td = 1052 and a = 0.303. (d) The 
dependence of the detrapping time Td on the number L of 
initially excited sites in log-log scale. The dashed straight 
line corresponds to a function oc L'*. 



to some diffusion process, and is characterized by some 
momentary diffusion rate D{t) such that 7712 — D{t)t. 
Then it follows that D{t) - t^^/s .^^^ finally D - 77^. 
Such a result has to be the outcome of the action of the 
nonlinear terms (which always contain products (3n). A 
diffusion rate is equal to an inverse characteristic time 
scale, and therefore we conjecture 



D = T- 



/3V 



(17) 



There are two ways of modifying D. We can either spread 
our initial excitation over some number of sites L, there- 
fore varying n. Alternatively we can fix the shape of the 
initial excitation, and vary (3. 

In order to test the validity of Eq. (fTT]) for a fixed value 
of nonlinearity we considered a single site excitation in 
the intermediate regime II for the KG model with total 
energy E = 0.4, so that 777,2 and P start to grow from 
the beginning (black curves in Figs. [8] (a) and (b) re- 
spectively). We also followed the time evolution of wave 
packets having as initial condition a homogeneous dis- 
tribution of the energy E — 0.4 among L neighboring 
sites. In particular, we considered initial distributions 
with 7t; = and p; = except for the central L sites 
whose initial momenta were set to ±yj2E/L, with the 
sign changing randomly from site to site. We performed 
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simulations with L ranging from 1 up to 41. The time 
evolution of m2 and P for some of these cases is shown 
in Figs. [H^a) and (b) respectively. In accordance to the 
results presented in Fig. [7] we observe that, distributing 
the energy of a single site excitation belonging to regime 
II over more sites results in a time dependence of m2 
and P similar to regime I, i. e. both quantities start to 
increase after some transient detrapping time r^. 

The behavior of the second moment m2 (t) can be mod- 
eled by a function of the form 

TO2(i) =M(i + rd)", (18) 

where M is a constant related to the value of the sec- 
ond moment of the initial distribution m2(0) = Mt^. 
Eq. ()18p gives a power law dependence of 7712 on t for 
t ^ Td and a slow time dependence of m2 for t <^ Td- 
Thus, it can be used to describe the behavior of m2 for 
L > 1. Fitting the numerical data obtained for different 
values of L by Eq. (see Fig.[5Kc) for such an example) 
we can determine the dependence of Td on L (Fig. [SKd)). 
Since L ^ from (fT7|) we conclude that Td ^ L'^- As 
we can see from Fig. [SJd) the numerically obtained re- 
sults are in good agreement with this assumption. 

To test the dependence of D on /?, we studied the weak 
nonlinearity regime I for the DNLS model with W — 4. 
We launched single site excitations for 10 realizations for 
(3 = 0.1 and /3 — 0.2. We estimated the detrapping 
times Td on logarithmic scale for each run, and averaged 
over each group of realizations. As a result we obtain 
(logioTd) = 5 for /? = 0.2, and (logig r^) = 6.9 for 
P = 0.1 (with (• • • ) denoting the mean value over the 
realizations), and their difference is then 1.9. According 
to (fT7)) . the difference should be 1.2 which is in relatively 
good agreement with the numerically estimated value. 

F. Numerical accuracy and roundoff errors 

We performed several tests in order to ensure that our 
results are not generated by inaccurate computations. 
First we varied the size of the system and found no de- 
pendence of the results on it. Therefore we exclude finite 
size effects. 

Second we varied the time steps of the symplectic in- 
tegration schemes by orders of magnitudes. Again we 
found no visible change in the detrapping times, or in 
the spreading characteristics. We also used difl^erent inte- 
gration schemes, and even nonsymplectic ones (8th order 
Runge-Kutta). No changes were obtained either. There- 
fore we exclude effects due to discretization of time. 

Finally we studied the influence of computational 
roundoff errors. The above observation, that the vari- 
ation of time steps does not change the key results, im- 
plicitly tells that roundoff errors can be excluded as well. 
Indeed, changing the time steps, we change the number 
of operations to be performed on a given interval of in- 
tegration. Therefore we change the number of roundoff 
operations. 



In addition, we decided to perform further tests with 
respect to the roundoff error issue. These tests are in- 
spired by the following consideration. Floating point 
numbers are characterized by the number of digits a after 
the comma which are kept during computations. All pre- 
sented data were obtained with double precision, where 
a = 16. The detrapping and spreading can be only due 
to the cubic nonlinear terms in the equations of motion. 
These terms are added to linear terms, when calculating 
the rhs of ^ and ([7]). Therefore, when for example in the 
case of the DNLS model sup; \tpi\'^ < 10~^, the nonhnear 
terms become of the order of the roundoff error of the 
linear terms. For all of our simulations, the amplitudes 
in the packet are of the order of 10~^ or larger. Therefore 
the roundoff is affecting only the amplitudes far in the 
exponential tails. We changed the calculation to single 
precision, for which a = 8, but we did not observe any 
qualitative difference in our results. For single precision 
the nonlinear terms will be affected by roundoff errors 
when sup; < 10"'', which is still realized only in 
the exponential tails. We note, that the times at which 
the roundoff errors affect the packet modes correspond 
to t - 10^° for a = 16 and t - lO^" for a = 8 which are 
obviously not accessible with our computation schemes. 

Therefore we implemented a brute force roundoff 
scheme: after each time step of integration we take the 
distributions and perform a roundoff at a prescribed digit 
a = l,2,3,4,.... We expect therefore to reduce the time 
at which roundoff errors will become visible, in order to 
observe that effect within the time window accessible by 
our computations. Indeed, we find that strong fluctu- 
ations in the conserved quantities set in at a time tr 
which decreases with decreasing a. In particular for the 
DNLS we find U « 10^10^10^ for a = 1,2,3, and for 
the KG model we find U « 10^ 10^ 10* for a = 1, 2, 3. 
When monitoring the second moment and the partici- 
pation number, we also find strong deviations from the 
above results at times t > tr- For a > 4 we do not ob- 
serve any significant change in the data. Therefore we 
conclude, that the roundoff errors with double (or even 
single) precision are not affecting our results. 



IV. SPREADING MECHANISMS 

We can think of two possible mechanisms of wave 
packet spreading. A NM with index fi in a, layer of width 
p in the cold exterior, which borders the packet, is ei- 
ther incoherently heated by the packet, or resonantly ex- 
cited by some particular NM from a layer with width p 
inside the packet. Heating here implies a (sub)diffusive 
spreading of energy. Note that the numerical results yield 
subdiffusion, supporting the nonballistic diffusive heating 
mechanism. 

For heating to work, the packet modes <j>i,{t) should 
contain a part (j)'^ {t) , having a continuous frequency spec- 
trum (similar to a white noise), in addition to a regular 



part (t)l,(t) of pure point frequency spectrum: 



(19) 



Therefore at least some NMs of the packet should evolve 
chaotically in time. The more the packet spreads, the 
less the mode amplitudes in the packet become. There- 
fore its dynamics should become more and more regular, 
implying limt_oo (j>l[t) / (plit) 0. 



A. Are all packet modes chaotic? 

In Ref. [8| it was assumed that all NMs in the packet 
are chaotic, and their phases can be assumed to be ran- 
dom at all times. At variance to the above expectation, it 
follows that 0^(i) = 0, or at least the ratio 4>1{t) / (f)l,{t) is 
constant on average. Consequently |(^^(t)| ~ n^/^ where 
n is the average norm density in the packet. 

According to ^ the heating of the exterior mode 
should evolve as ~ A^^^ + f3n-^/^f{t) where 

{f{t)f{t')) = ^{t — t') ensures that f{t) has a continu- 
ous frequency spectrum. Then the exterior NM increases 
its norm according to \4>^\'^ ^ 0^TV't. The momentary 
diffusion rate of the packet is given by the inverse time 
T it needs to heat the exterior mode up to the packet 
level: D = 1/T ^ 0^1-1^ . The diffusion equation m2 ^ Dt 
yields ra2 ~ pt^^'^. We tested the above conclusions by 
enforcing decoherence of NM phases. Each 100 time units 
on average 50% of the NMs were randomly chosen, and 
their phases were shifted by tt (DNLS). For the KG case 
we changed the signs of the corresponding NM momenta. 
We obtain m2 ~ t^^"^ (see FigEl). Therefore, when the 
NMs dephase completely, the exponent a = 1/2, contra- 
dicting numerical observations without dephasing. Thus, 
not all NMs in the packet are chaotic, and dephasing is 
at best a partial outcome. 



B. Mode- mode resonances inside the packet 

Chaos is a combined result of resonances and noninte- 
grability. Let us estimate the number of resonant modes 
in the packet for the DNLS model. Excluding secular 
interactions, the amplitude of a NM with = Ui, is 

modified by a triplet of other modes fl = (/ii, /12, /is) in 
first order in as (|4]) 



I — /^V^Mi '^Ma a ' 



dX 



1^1,(^2,(^3 



, (20) 



where dX 



X^ + X, 



A„ 



The perturba- 



tion approach breaks down, and resonances set in, when 
< {(ji'^J'^l- Since all considered NMs belong to the 
packet, we assume their norms to be equal to n for what 
follows. If three of the four mode indices are identical, 
one is left with interacting NM pairs. A statistical anal- 
ysis of the probability of resonant interaction was per- 
formed in Ref. JTj. For small values of n (i.e. when the 




FIG. 9: (color online) Statistical properties of NMs of the 
DNLS model. Prob ability densities W{Ru.Po) of NMs being 
resonant (see sect ion HV Bl for details) . Disorder strength W = 
4, 7, 10 (from top to bottom). 



packet has spread over many NMs) the main contribution 
to resonances are due to rare multipeak modes [ll| , with 
peak distances being larger than the localization volume. 
If two or none of the four mode indices are identical, one 
is left with triplets and quadruplets of interacting NMs 
respectively. In both cases the resonance conditions can 
be met at arbitrarily small values of n for NMs from one 
localization volume. 

We perform a statistical numerical analysis for the 
quadruplet case. For a given NM v we obtain R^.p„ — 
ininp R,y p. Collecting Ri,^p(, for many v and many dis- 
order realizations, we find the probability density dis- 
tribution W{R^jg) (Fig. ini). The main result is that 
y^{Ru,iio -> 0) ^ C{W) ^ 0. For the cases studied, the 
constant C drops with increasing disorder strength W . 
Similar results are found if pairs of resonant NMs |Ti| 
are analyzed, with the only difference that the constant 
C is reduced e.g. by a factor of 30 for W = 4. 

The probability V for a mode, which is excited to a 
norm n (the average norm density in the packet), to be 
resonant with at least one triplet of other modes at a 
given value of the interaction parameter /3 is given by 



V = 



W{x)<lx , 



with X denoting R^^p.^ ■ For /3n ^ 1 it follows 



(21) 



(22) 



Therefore the probability for a mode in the packet to be 
resonant is proportional to C(3n. On average the num- 
ber of resonant modes in the packet is constant, pro- 
portional to C/3, and their fraction within the packet is 
^ C(3n. Since packet mode amplitudes fluctuate in gen- 
eral, averaging is meant both over the packet, and over 
suitably long time windows (yet short compared to the 
momentary inverse packet growth rate). We conclude, 
that the continuous frequency part of the dynamics of a 
packet mode is scaled down by C'(3n, compared to the 
case when all NMs would be chaotic. It follows that 
4>'^{t)/(l)l,{t) ~ C/3n. As expected initially, the chaotic 



10 



part in the dynamics of packet modes becomes the weaker 
the more the packet spreads, and the packet dynamics be- 
comes more and more regular in the hmit of large times. 
Therefore the chaotic component <j>l{t) <C <P^{t) is a small 
parameter. 

Expanding the term ipi, to first order in (j>'^{t), the 
heating of the exterior mode should evolve according to 
« + It follows _ c^p^nH, 

and the rate D = l/T - C^/J^n^ (cf. the prediction (fT7| V 
The diffusion equation 7712 ~ Dt yields 

m2 - C72/3/34/3t" , a = 1/3 . (23) 

The predicted exponent a — 1/3 is close to the numeri- 
cally observed one, as we discussed in section ITlI Dl 

C. Resonant spreading? 

Finally we consider the process of resonant excitation 
of an exterior mode by a mode from the packet. The 
number of packet modes in a layer of the width of the 
localization volume at the edge, which are resonant with 
a cold exterior mode, will be proportional to /?n. After 
long enough spreading /3rt 1. On average there will 
be no mode inside the packet, which could efficiently res- 
onate with an exterior mode. Therefore, resonant growth 
can be excluded. 



V. SUMMARY AND DISCUSSION 

We studied the spreading of wave packets in disordered 
one-dimensional nonlinear chains. In particular we con- 
sidered two systems, namely the DNLS model (HJ and 
the quartic KG system ^ . The linear parts of these two 
models are equivalent in the sense that they correspond 
to the same eigenvalue problem ([3]). 

We predicted theoretically and verified numerically the 
existence of three different dynamical behaviors depend- 
ing on the relation of the nonlinear frequency shift S 
(which is proportional to the system's nonlinearity) with 
the average spacing AA of eigenfrequencies and the spec- 
trum width A (AA < A) of the linear system. The dy- 
namics for small nonlinearities {S < AA) is characterized 
by localization as a transient, with subsequent subdiffu- 
sion (regime I). For intermediate values of the nonlinear- 
ity AA < 5 < A, and the wave packets exhibit immediate 
subdiffusion (regime II). In this case, the second moment 
7712 and the participation number P increase in time fol- 
lowing the power laws m2 ~ t", P ~ i"/^. Assuming 
that the spreading is due to an incoherent excitation of 
the cold exterior, induced by the chaotic behavior of the 
wave packet, we predicted a = 1/3. Finally, for even 
higher nonlinearities ((5 > A) a large part of the wave 
packet is selftrapped, while the rest subdiffuses (regime 
III). In this case P remains practically constant, while 
7712 ~ t"- The overall picture is schematically presented 
in Fig. [U both for the DNLS and the KG model. 



The compactness index C = P^/m2, which measures 
the sparseness of wave packets, exhibits different be- 
haviors for the three dynamical regimes. In particular, 
the behavior of C for wave packets in regime II imply 
that these wave packets spread but do not become more 
sparse. 

For large values of the disorder strength W and/or 
strong nonlinearity the intermediate regime II effectively 
disappears, and the evolution will start either in regime I, 
or in regime III. In regime I the detrapping times increase 
with further increase of W. In regime III the fraction of 
the wave packet which spreads decreases with increas- 
ing nonlinearity. Therefore, large values of W and/or 
nonlinearity will not allow for an observation of the de- 
struction of Anderson localization on time scales which 
are bounded from above by practical computational lim- 
itations. 

The subdiffusive spreading is universal, i. e. the expo- 
nent a is independent of the nonlinearity's strength (/3 
for the DNLS model and energy E for the KG one) and 
W, which are only affecting the prefactor in Ex- 
cluding selftrapping, any nonzero nonlinearity strength 
P will completely delocalize the wave packet and destroy 
Anderson localization. The exponent a is determined 
solely by the degree of nonlinearity, which defines the 
type of overlap integral to be considered in ^U\i . and by 
the stiffness of the spectrum {A,y}. Our numerical com- 
putations confirmed the prediction a = 1/3 in the case 
of single site and of nonlocal homogeneous excitation. In 
the case of single mode excitations the three different 
regimes were also detected. The numerically computed 
exponents a get slightly larger values than 1/3, exhibit- 
ing also a small dependence on the strength of nonlin- 
earity. This discrepancy between the two cases in not 
clearly understood. 

We studied the statistics of detrapping times for 
regime I. We provided numerical evidences for the valid- 
ity of the conjectured dependence of Td on the nonlin- 
earity strength and on the average norm density of the 
excited NMs given in Eq. ^T7\ . It is worth mentioning 
that, distributing the energy of a single site excitation 
belonging to regime II over more sites results in a time 
dependence of 1112 and P similar to regime I. In addition, 
considering as initial condition the profile of a single site 
excitation in regime II at some latter time td, we observe 
a dynamical evolution of the type of regime I where the 
detrapping time is Td ^ td- 

The spreading of the wave packet is due to weak but 
nonzero chaotic dynamics inside the packet. It is natural 
to expect such a dynamics, since the considered systems 
are nonintegrable. If instead an integrable system is con- 
sidered, Anderson localization will not be destroyed. In- 
deed, consider a Hamiltonian in NM representation using 
actions Ji, and angles 9^, as coordinates: 



(24) 
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We assume that the set of eigenfrequencies {X,^} and the 
overlap integrals /, 



1^1 ,1^2 ,1^3 ,1^4 



are identical with those de- 
scribing the DNLS model (H]) , ([5|) . The equations of mo- 
tion — —dHint/dOu and 9u — dTimt/dJu yield Ju = 
since the integrable Hamiltonian (I24p depends only on 
the actions. Therefore, any localized initial condition 
(e. g. J^(t = 0) oc S^,j^g) will stay localized, since ac- 
tions of modes which are at large distances will never get 
excited. Thus, the observed spreading of wave packets, 
which we studied in detail in the present work, is entirely 
due to the nonintegrability of the considered models, at 
variance to ([24]) . 

The more the wave packet spreads, the weaker the res- 
onances become. Corresponding structures (chaotic lay- 
ers) in phase space become thinner and thinner. Consider 
quantum many-body systems. Classical phase space 
structures which are finer than the action quantization 
induced grid become irrelevant. Therefore we may spec- 
ulate, that the wave packet will stop spreading for a quan- 
tum many-body system at some point for zero tempera- 
ture, but also for temperatures below some finite thresh- 
old. These expectations are very close to rigorous results 
for interacting fermions in disordered systems (2^ . 

In our study we considered initial conditions exciting 
NMs whose eigenvalues are located close to the center of 
the frequency band. Thus, the evolution of the system 
does not significantly depend on the sign of nonlinear- 
ity. In contrast, when one excites eigenstates with fre- 
quencies near the band edges, a rather weak nonlinearity 
might lead either to selftrapping or to the weak nonlin- 
ear regime depending on the sign of nonlinearity. Such 
examples were presented in ^] where NMs close to the 
edges of the band exhibit different dynamical behaviors, 
i. e. one becomes more localized as the nonlinearity was 
switched on, while the other tends to delocalize. If a 
spatially continuous system is considered, then a proper 
choice of the sign of nonlinearity prohibits selftrapping 
(so-called defocusing nonlinearity, corresponding to re- 
pulsive two-body interactions). For such a case, regime 
III ceases to exist, and localization is expected to be de- 
stroyed irrespectively of the strength of nonlinearity. 
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APPENDIX A: THE SABAa AND SBABa 
SYMPLECTIC INTEGRATORS 

In a family of symplectic integrators which involve 
only forward integration steps was proposed. These inte- 
grators were adapted for integrations of perturbed Hamil- 
tonians of the form 



where both A and B are integrable and e is a parameter. 
We briefly recall here their main properties focusing our 
attention on two particular members of the family of inte- 
grators presented in [l8| , namely the S ABA2 and SBAB2 
integrators. These integrators have already proved to be 
very efficient for the numerical study of astronomical , 
as well as accelerator models '21]. 

Consider a Hamiltonian system of N degrees of 
freedom having a Hamiltonian H{p,u), with p ~ 
{pi,...,Pn), u = (ui,...,UAr) where ui and pi, I = 
1, . . . , iV, are the generalized coordinates and momenta 
respectively. An orbit of this system is defined by a vec- 
tor x{t) = {xi{t), . . .,X2N{t)), with Xi = pu Xi+N = Ui, 

I — 1, . . . , N. This orbit is a solution of Hamilton's equa- 
tions of motion: 



dpi 
dt 



dH 
dui 



dui 
dt 



dH 
dpi 



.,7V, (A2) 



where t is the independent variable, namely the time. 
Defining the Poisson bracket of functions /(p, u), g{p, u) 
by: 



^^ydp^du^ dui dpi J' 
the Hamilton's equations of motion take the form: 



(A3) 



dx 

— = {H,x} = Lhx, 



(A4) 



where Lh is the differential operator defined by L^/ = 
{x, /}. The solution of Eq. (|A4[) . for initial conditions 
a;(0) — xq, is formally written as: 



f 

z — ^ 77 1 



Xo- 



(A5) 



ri>0 



A symplectic scheme for integrating (jA4p from time t 
to time t + T consists of approximating, in a symplectic 
way, the operator e'^^" — e'^^^^"'"^^^-' by an integrator 
of j steps involving products of e'^*'^'^-* and e'*'^^''^, i — 
l,2,...,j, which are exact integrations over times 
and diT of the integrable Hamiltonians A and B. The 
constants Ci, di, are chosen so that to increase the order 
of the remainder of this approximation. 

For the SABA2 integrator we get: 

SABA2 = e^'^LA^d.rL.B^c^TLA^d.rL.B^c.rLA^ (^g) 

with ci = i(^l-^^, C2 = di = i, while the 
SBAB2 integrator is given by 



SBAB^ 



(A7) 



H = A + eB, 



(Al) 



with C2 = di = ^, ^2 = |- Using these integrators 
we are actually approximating the dynamical behavior 
of the real Hamiltonian A + ei? by a Hamiltonian H = 
A + eB + 0{T*e + r^e^), i. e. we introduce an error term 
of the order r^e -t- r^e^. 
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The accuracy of the SABA2 (or SBAB2) integrator can 
be improved when the term C — {{A, leads to 

an integrable system, as in the common situation of A 
being quadratic in momenta p and B depending only on 
positions u. In this case, two corrector terms of small 
backward steps can be added to the integrator SABA2 



SABA2C 



^2 9 r ^ 



(SABA2)e' 



(A8) 



A similar expression is valid also for SBAB2 . The value of 
g was chosen in order to eliminate the r^e^ dependence of 



the remainder which becomes of order 0{t e - 



^). In 



particular we have g = (2 — \/3)/2 for SABA2 and g = j2 
for SBAB2. We note that the SABA2 and SBAB2 inte- 
grators involve only forward steps which increases their 
numerical stability, while, the addition of the corrector 
results to better accuracy of the schemes, introducing si- 
multaneously a small backward step. 



with: 



N 



A = \ ^ 



1=1 

N 



B = 



^ 2 

1, 



1 4 

-uf 
4 ' 



uif 



(A9) 



where N is the number of anharmonic oscillators. The 
operators e^^^, e'^^^, e^^'^ , which propagate the set of 
initial conditions {ui,pi) at time t, to their final values 
{u[,p[) at time t + t, I ^ 1,2, . . . , N are: 



1. Integration of the KG lattice 

Hamiltonian ([6]) is suitable for the implementation of 
the SABA2C integration scheme since it attains the form 



Ui = Pit 
P'l = Pi 



ui 



(AlO) 




-ui (e/ + ) + — + Ui+i - 2ui) 



T+Pl 



(All) 



p'l 



Pi = 



Pn 



ui 
2 



W 



ei + 3uf 



-Ml (ei +ul) + — {u2 - 2ui) 



W 

■A- 

I w 



U2 (e2 +uf) ~ — (W3 + Ml - 2li2) 



T +Pl 



+ 



+ 



ui-1 (e;-i + uf_j) - — (ii;_2 +ui~ 2ii/_i) 



-ui {li + ) + — (u/_i + ui+i - 2ui) 



W 



T + pi, for ? = 2,3, ... ,7V- 1 



W 
2 



+ I — -f ?jv + 



-UN (cat + U^) + — (mat-i - 2uAr) 



T ^Pn 



(A12) 



since 



N 



1=1 



ui (ei + uf) - — (u/_i + U/+1 - 2u/) 



(A13) 



and uq = ujv+i — 0. 



2. Integration of the DNLS lattice 



nian ([T]) as 



We use the SBAB2 integrator scheme to integrate the 
equations of motion ^ , by splitting the DNLS Hamilto- 



N 



^ = -^i'>Pi+ii^i +^i+i'^i)' 



1=1 



N 



1=1 



(A14) 



13 



with being the number of lattice sites. The action of 
the operator e^^^ on tpi, I = 1,2, . . . , N at time t leads 
to the computation of V; at time t + t, and includes 
three steps: a) the transformation of the wavefunction 
from the real {ipi) to the Fourier {(pq) space, through a 
Fast Fourier transform (FFT), b) a rotation of (pq, and 
c) the inverse FFT of the wavefunction Pq evaluated at 
the previous step, i. e. 



-J^^ : < 



2Txiq(m-l)/N 



^/^ = ^^g2icos(27r(g-l)/JV)r 



(A15) 



1 

N ^1 



E;=i^;e-2-K,-i)/iv 



On the other hand, the action of e^^^ on reduces to 
a simple rotation in real space, namely 



Note that for the DNLS model we do not apply the two 
corrector steps since the term C = {{A, B), B) does not 
lead to an easily solvable system. 



(A16) 
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